/*
**  File: AquaFuncts_pdb.c
**
**  This file is part of the AQUA program suite, which is being
**  developed in the NMR Spectroscopy Department, Bijvoet Center for 
**  Biomolecular Research, Utrecht. 
**
**
**  AUTHOR(S) : Ton Rullmann
**  VERSION   : 1.0
**  DATE      : Nov 5, 1996
**
**  This file contains routines for handling PDB files.
**
**  Most header definitions are present in AquaFuncts_io.c
*/

/****************  function prototypes  ****************/

#include "AquaFuncts_pdb.h"

static int ReadHeader_pdb( FILE *fp, char *header, const size_t size );
static Status ReadAtomData_pdb( FILE *fp, char *atnam, char *resnam, char *chn, 
                                char *ins, int *atnum, int *resnum, 
			        const int size, const int count );
static Status ReadCoordData_pdb( FILE *fp, float *x, float *y, float *z );
static Status ReadElementData_pdb( FILE *fp, char *el );
static void GuessElement_PDB( char *a, char *r, char *el );
static void FormatElement_PDB( char *el );

/* BEGIN_EXTH */

/****************  ReadFil_cv_pdb  ****************/
/*
**  Read a PDB file in the ReadFil_cv way.
*/

Status ReadFil_cv_pdb( FILE *fp, const Filtyp type, const int atom_count,
                       int *fil_struc_num, int *nerr )
{
    int      count = 0;
    char     *a, *r, *chn, *ins, *el;
    int      *ia, *ir;
    float    *px, *py, *pz;
    Status   stat;

/* allocate memory */

    Alloc_cv( atom_count );
        
/* read the file header - for the first structure */
/* or the ENDMDL/MODEL records - for subsequent structures */
/* NOTE: ReadHeader_pdb returns the MODEL number read..... */


    *fil_struc_num = ReadHeader_pdb( fp, header, size_HEADER );
    
/* read the predetermined number of atoms */

    a = atonam;
    r = resnam;
    chn = chain;
    ins = insert;
    ia = atono;
    ir = resno;
    px = x;
    py = y;
    pz = z;
    el = atoel;

    while ( ++count <= atom_count )
    {
        stat = ReadAtomData_pdb( fp, a, r, chn, ins++, 
	                         ia++, ir++, size_ID, count );
        if ( stat != status_ok )
	{
	    if ( stat == eof_err && count == 1 )
		return( end_of_data );
	    else
                goto err;
	}
        stat = ReadCoordData_pdb( fp, px++, py++, pz++ );
        if ( stat != status_ok )
            goto err;

        stat = ReadElementData_pdb( fp, el );
        if ( stat != status_ok )
            goto err;
        
        if ( OptionTrue( opt_element_guess_pdb ) ) {
            GuessElement_PDB( a, r, el );            
        }
                
        FormatElement_PDB( el );
                
        a   += size_ID;
        r   += size_ID;
        el  += size_ID;
	chn += size_ID;

        if ( FlushLine( fp ) < 0 )
	{
	    stat = eof_err;
            goto err;
	}
    }
    return( status_ok );

err:
    *nerr = count-1;
    puts( " * Error or unexpected end of file" );
    return( stat );
}

/****************  AtomCount_pdb  ****************/
/*
**  Count the number of atoms per structure in a PDB file.
*/

int AtomCount_pdb( FILE *fp )
{
    Card     line;
    int      count = 0;
    Boolean  model = false;
    
    while ( ReadLine( line, size_LINE, fp ) )
    {
        if ( EQn( line, "ATOM", 4 ) || EQn( line, "HETATM", 6 ) )
	{
	    count++;
	    continue;
	}
	if ( EQn( line, "ENDMDL", 6 ) || EQn( line, "END", 3 ) ||
	     EQn( line, "CONECT", 6 ) || EQn( line, "MASTER", 6 ) )
	    break;
	if ( EQn( line, "MODEL", 5 ) )
	{
	    if ( model )
		break;
	    else
		model = true;
	}
    }
    return( count );
}

/****************  Write_pdb  ****************/
/*
**  Write a PDB file, either a single model: the model with number <model>, 
**  or all models: if <model> equals 0.
**  In case of a single model, MODEL/ENDMDL records are written if
**  <wrmodel> is true. For multiple models, MODEL/ENDMDL records are always
**  written.
*/

Boolean Write_pdb( FILE *fp, Strucset *set, const int model, Boolean wrmodel )
{
    Struc   *strucptr;
    Coord   *coordptr;
    Mol     *molptr = set->molecule;
    Atom    *atomptr;
    int     count, num = 0;
    ID      temp_id;
        
    if ( NOTEQn( set->ident, "REMARK", 6 ) )
	fprintf( fp, "REMARK " );
    fprintf( fp, "%s\n", set->ident );
    fprintf( fp, "REMARK %s\n", set->derived );

    if ( set->struc_count > 1 && !model )
	wrmodel = true;
    strucptr = set->fst_struc;
    while ( strucptr )
    {
	if ( !model || strucptr->struc_num == model )
	{
	    if ( wrmodel )
		fprintf( fp, "MODEL     %4i\n", strucptr->struc_num );
	    count = set->atom_count;
	    coordptr = strucptr->coords;
	    atomptr = molptr->atoms;
	    while ( count-- )
	    {
		if ( num == 99999 )
		    num = 0;
		fprintf( fp, "ATOM  %5i ", ++num );
		if ( LenChars( atomptr->atom_ext, size_ID ) == 4 ||
		     isdigit( atomptr->atom_ext[0] ) )
		{
		    FPrintChars( fp, atomptr->atom_ext, 4 );
		}
		else
		{
		    fputc( ' ', fp );
		    FPrintChars( fp, atomptr->atom_ext, 3 );
		}
		fputc( ' ', fp );
		FPrintChars( fp, atomptr->res_nam, 3 );
		fprintf( fp, " %c%4i%c   %8.3f%8.3f%8.3f%6.2f%6.2f          ", 
		    atomptr->chain[0], atomptr->res_num_orig, atomptr->res_insert, 
		    coordptr->x, coordptr->y, coordptr->z, 1.0, 0.0);
                
		CopyChars(temp_id, atomptr->atom_element_symbol, 2);
                temp_id[1] = toupper(temp_id[1]);
                temp_id[2] = '\0';
                if ( temp_id[1] == ' ' ) {
                    temp_id[1] = temp_id[0];
                    temp_id[0] = ' ';
                }
                    
		fprintf( fp, "%2s\n", temp_id);
                
		coordptr++;
		atomptr++;
	    }
	    if ( wrmodel )
		fprintf( fp, "ENDMDL\n" );
	}
        strucptr = strucptr->next_struc;
    }

    return( true );
}


/* END_EXTH ----------------------- only private functions follow */

/****************  ReadHeader_pdb  ****************/
/*
**  Read header information from a PDB file
**  Format assumed: 
**  - text header copied from lines beginning with: HEADER, COMPND, AUTHOR, 
**    REVDAT or REMARK
**    (only the first MAXHEADER of such lines are stored).
**  - the function returns after reading a SCALE3 or MODEL record.
**  - if a MODEL record is encountered, the function returns the model
**    number, otherwise 0.
**  - if an ATOM or HETATM record is encountered before a SCALE3 or
**    MODEL record, the file is rewound.
**
**  After returning the file is positioned at the beginning of a new line.
*/

static int ReadHeader_pdb( FILE *fp, char *header, const size_t size )
{
    Card  line;
    char  *h = header;
    int   numh = 0, model = 1;

    while ( ReadLine( line, size_LINE, fp ) )
    {
        if ( EQn( line, "HEADER", 6 ) || EQn( line, "COMPND", 6 ) || 
             EQn( line, "AUTHOR", 6 ) || EQn( line, "REVDAT", 6 ) ||
	     EQn( line, "REMARK", 6 ) )
        {
            if ( numh++ < MAXHEADER )
            {
                strncpy( h, line, size );
                h += size;
            }
        }
        else if ( EQn( line, "SCALE3", 6 ) )
            break;
        else if ( EQn( line, "MODEL", 5 ) )
	{
	    sscanf( &line[6], "%i", &model );
            if ( numh++ < MAXHEADER )
            {
                strncpy( h, line, size );
                h += size;
            }
            break;
	}
        else if ( EQn( line, "ATOM", 4 ) || EQn( line, "HETATM", 6 ) )
        {
	    puts( "Warning - rewinding the input file before reading next structure" );
            rewind( fp );
            break;
        }
    }
    return( model );
}

/****************  ReadAtomData_pdb  ****************/
/*
**  Read atom information from one PDB ATOM or HETATM record
**  Format assumed: (read up to and including column 27)
**  123456789012345678901234567890123456789012345678901234567890123456789012345678901
**           1         2      X  3         4         5         6         7
**  ATOM      1  H1  MET     1      27.340  26.800   6.430  1.00  0.00
**  ATOM    496 1HH2 ARG A  51      30.100  21.310  13.800  1.00  0.00    
**  HETATM  748 ZN    ZN    47       6.128  -4.825   1.297  1.00  0.00          ZN
**
**  FORMAT(6A1,I5,1X,A4,A1,A3,1X,A1,I4,A1,3X,3F8.3)
**
**  The routine finds the first ATOM or HETATM record and reads it
**  (finding is necessary because ReadHeader_pdb does not always position
**  the file at an ATOM/HETATM record, and also to skip intervening records
**  such as SIGATM, TER, MODEL and ENDMDL).
**  CONECT, MASTER and END records, and premature end-of-file 
**  generate an error condition.
**
**  After reading the file is **not** flushed to the next line.
*/

static Status ReadAtomData_pdb( FILE *fp, char *atnam, char *resnam, char *chn, 
                                char *ins, int *atnum, int *resnum, 
			        const int size, const int count )
{
    char  loc, dumm1, dumm2;
    char  id[7];
    
    while ( true )
    {
        if ( fscanf( fp, "%6s", id ) == EOF )
	    return( eof_err );
				     /* read id; stop if eof */
        if ( EQn( id, "ATOM", 4 ) || EQn( id, "HETATM", 6 ) )
            break;
	                             /* only ATOM and HETATM records are processed */
	if ( EQn( id, "END", 3 ) ||
             EQn( id, "CONECT", 6 ) || EQn( id, "MASTER", 6 ) )
            return( eof_err );
	                             /* stop at CONECT, MASTER, END */
        if ( FlushLine( fp ) < 0 )
            return( eof_err );       /* stop if no more lines */
    }

    if ( fscanf( fp, "%5i%1c%4c%1c%3c%1c%1c%4i%1c", 
                 atnum, &dumm1, atnam, &loc, resnam, &dumm2, chn, resnum, ins ) 
         < 7 )
    {
        puts( " * Error in reading atom data from PDB ATOM/HETATM record" );
        return( read_err );
    }
    
    if ( loc != ' ' )
    {
	puts( " * Don't know how to handle non-blank <alternate location> field" );
	printf( " * PDB ATOM/HETATM record, atom number %i\n", *atnum );
/*	exit( 1 );  */
    }

/* remove initial space(s) - if present - from atom and residue names */

    RemLeadSpc( atnam, 4 );
    RemLeadSpc( resnam, 3 );
/*
    if ( *atnam == ' ' )
    {
	char *a;
	int  i = 0;
	a = atnam + 1;
	while ( ++i < 4 )
	    *atnam++ = *a++;
	*atnam = ' ';
    }
*/
    return( status_ok );
}

/****************  ReadCoordData_pdb  ****************/
/*
**  Read coordinate information from one PDB ATOM record
**  Format assumed: free format floats
**
**  After reading the file is **not** flushed to the next line.
**  Returns status flag.
**  Read atom information from one PDB ATOM or HETATM record
**  Format assumed: (read up to and including column 54)
**  123456789012345678901234567890123456789012345678901234567890123456789012345678901
**           1         2       B 3         4         5   E     6         7
**  ATOM      1  H1  MET     1      27.340  26.800   6.430  1.00  0.00
*/

static Status ReadCoordData_pdb( FILE *fp, float *x, float *y, float *z )
{

    /* Updated to make sure that this part ends with reading the char in column 54. */
    
    if ( fscanf( fp, "%*3c%8f%8f%8f", x, y, z ) < 3 )
    {
        puts( " * Error in reading atom data from PDB ATOM/HETATM record" );
        return( read_err );
    }
            
    return( status_ok ); 
}

/****************  ReadElementData_pdb  ****************
**  Returns status flag.
**  Format assumed: nothing needs to be there, if there is something, it will be stored
**  123456789012345678901234567890123456789012345678901234567890123456789012345678901
**           1         2         3         4         5    B    6         7       E
**  HETATM  748 ZN    ZN    47       6.128  -4.825   1.297  1.00  0.00          ZN
**  HETATM  748 ZN    ZN    47       6.128  -4.825   1.297  1.00  0.00            
**  HETATM  748 ZN    ZN    47       6.128  -4.825   1.297
*/

static Status ReadElementData_pdb( FILE *fp, char *el )
{

    char dummy[25];
    strcpy(el,"  ");
    
    fgets( dummy, 23, fp );
    if ( dummy[strlen(dummy)-1] == '\n' ) {
        ungetc( dummy[ strlen(dummy)-1 ], fp);
    } else {
        fgets( dummy,  5, fp );
        if ( dummy[strlen(dummy)-1] == '\n' ) {
            ungetc( dummy[ strlen(dummy)-1 ], fp);
        }
        
        if ( strlen(dummy) >= 2 ) {
            el[0] = dummy[0];               
            el[1] = dummy[1];               
        }

    }    
    return( status_ok ); 
}

/* This routine GuessElement_PDB should only be invoked when the atom names are
corrected to IUPAC and even then no garantees are given that it will always be
producing the correct output */

static void GuessElement_PDB( char *a, char *r, char *el)
{
    
    
    /* Return if el is not empty */
    if ( ! (( el[0] == ' ' ) && ( el[1] == ' ' )) ) {
        return;        
    }
    
    /* Make the element the same as the first character of the atom name */
    el[0] = a[0];
    
    
    /* Special cases:     ions
    if residue name and atom name share the first two characters */
    
    if ( ( a[0] == r[0] ) && ( a[1] == r[1] ) ) {
        el[0] = a[0];
        el[1] = a[1];   
    }
    
    /* Other cases here */
    /* When this becomes important we'll make a lib */

    /* Add the second char for some metals */
    if ( ( strncmp("CCH", r, 3) == 0 ) && ( strncmp("FE", a, 2) == 0 ) ) {
        el[1] = a[1];   
    }
    
}        


static void FormatElement_PDB( char *el)
{
    /* atom element symbol formatting*/

    RemLeadSpc( el, 2 );
    el[0] = toupper( el[0] );
    el[1] = tolower( el[1] );
}    
